Environment Setting

# Import required packages
library(tidytransit)
library(tidyverse)
library(tmap)
library(ggplot2)
library(gtfsrouter)
library(here)
library(units)
library(sf)
library(leaflet)
library(tidycensus)
library(plotly)
library(igraph)
library(tidygraph)
library(dodgr)
source("https://raw.githubusercontent.com/BonwooKoo/gtfs_to_igraph/master/gtfs_to_igraph.R")

wd <- file.path(Sys.getenv('setwd'),"work/working/School/UA_2022/external/Lab/module_2")
setwd(eval(wd))

What is the General Transit Feed Specification (GTFS)?

Before GTFS, there wasn’t (to the best of my knowledge, at least) a standardized format for transit timetables other associated information. Users of transit data across multiple transit agencies had to deal with different data formats. With GTFS, feeds from different agencies became standardized (although not perfect), making applications of the data easy.

Useful Links:

• tidytransit (i.e., R package) Vignettes
• Google GTFS
• Open Mobility Data - Archive of GTFS feeds from around the world

Let’s download some GTFS feed provided by Metropolitan Atlanta Rapid Transit Authority (MARTA).

## Directly from URL into R
# atl <- tidytransit::read_gtfs("http://www.itsmarta.com/google_transit_feed/google_transit.zip")

atl <- read_gtfs('MARTA_GTFS_Latest_Feed.zip')

# Take a look
summary(atl)

Note that the code above downloaded the GTFS feed from an URL directly into your R environment. If you want to get a stable copy for your self, you can also save the GTFS feed into your hard drive first and then read it into R. See below.

# NOT RUN: From URL to hard drive, then into R
download_path <- file.path(getwd(), "GTFS_MARTA.zip")
download.file("http://www.itsmarta.com/google_transit_feed/google_transit.zip", destfile = download_path)
atl <- tidytransit::read_gtfs(download_path) 

GTFS feed contains many relational tables about transit service schedules, trips, stops, and routes.

Understand GTFS

GTFS consists of multiple tables and comes in a zip file as a single package. Transit is a complex system that contains multiple components (e.g., routes, stops, service schedules) working together. The table below shows a brief description of what each data.frame contains. This table is taken from Google.

Description of tables in GTFS feed Table name Defines ———————- ———————
agency Transit agencies with service represented in this dataset.
stops Stops where vehicles pick up or drop off riders. Also defines stations and station entrances. routes Transit routes. A route is a group of trips that are displayed to riders as a single service. trips Trips for each route. A trip is a sequence of two or more stops that occur during a specific time period. stop_times Times that a vehicle arrives at and departs from stops for each trip. calendar Service dates specified using a weekly schedule with start and end dates. This file is required unless all dates of service are defined in calendar_dates.txt. calendar_dates Exceptions for the services defined in the calendar.txt. If calendar.txt is omitted, then calendar_dates.txt is required and must contain all dates of service. fare_attributes Fare information for a transit agency’s routes. fare_rules Rules to apply fares for itineraries. shapes Rules for mapping vehicle travel paths, sometimes referred to as route alignments. frequencies Headway (time between trips) for headway-based service or a compressed representation of fixed-schedule service. transfer Rules for making connections at transfer points between routes. pathways Pathways linking together locations within stations. levels Levels within stations. feed_into Dataset metadata, including publisher, version, and expiration information. translations Translated information of a transit agency. attributions Specifies the attributions that are applied to the dataset. ———————- ———————

These tables are relational table that are connected through a system of join keys. The schematic below shows which tables are linked to which tables, through which join keys. Understanding this structure is essential in using GTFS.

IMAGE SOURCE: http://tidytransit.r-transit.org/articles/introduction.html

typeof(atl)

## [1] "list"

names(atl)

## [1] "agency"         "calendar"       "calendar_dates" "routes"        
## [5] "shapes"         "stop_times"     "stops"          "trips"         
## [9] "."

# head(atl$calendar)

GTFS into geospatial format

The function gtfs_as_sf converts ‘shapes’ and ‘stops’ tables in GTFS data into sf objects.

atlsf <- tidytransit::gtfs_as_sf(atl, crs = 4326)
head(atlsf)

## $agency
## # A tibble: 1 × 7
##   agency_id agency_name                  agenc…¹ agenc…² agenc…³ agenc…⁴ agenc…⁵
##   <chr>     <chr>                        <chr>   <chr>   <chr>   <chr>   <chr>  
## 1 MARTA     Metropolitan Atlanta Rapid … http:/… Americ… en      (404)8… custse…
## # … with abbreviated variable names ¹​agency_url, ²​agency_timezone,
## #   ³​agency_lang, ⁴​agency_phone, ⁵​agency_email
## 
## $calendar
## # A tibble: 4 × 10
##   service_id monday tuesday wednesday thursday friday saturday sunday start_date
##   <chr>       <int>   <int>     <int>    <int>  <int>    <int>  <int> <date>    
## 1 2               0       0         0        0      0        0      0 2021-08-14
## 2 3               0       0         0        0      0        1      0 2021-08-14
## 3 4               0       0         0        0      0        0      1 2021-08-14
## 4 5               1       1         1        1      1        0      0 2021-08-14
## # … with 1 more variable: end_date <date>
## # ℹ Use `colnames()` to see all variable names
## 
## $calendar_dates
## # A tibble: 3 × 3
##   service_id date       exception_type
##   <chr>      <date>              <int>
## 1 5          2021-09-06              2
## 2 5          2021-11-25              2
## 3 5          2021-11-26              2
## 
## $routes
## # A tibble: 118 × 8
##    route_id route_short_name route_lon…¹ route…² route…³ route…⁴ route…⁵ route…⁶
##    <chr>    <chr>            <chr>       <chr>     <int> <chr>   <chr>   <chr>  
##  1 15664    1                Marietta B… ""            3 https:… FF00FF  ""     
##  2 15665    2                Ponce de L… ""            3 https:… 008000  ""     
##  3 15666    3                Martin Lut… ""            3 https:… FF8000  ""     
##  4 15667    4                Moreland A… ""            3 https:… FF00FF  ""     
##  5 15668    5                Piedmont R… ""            3 https:… 00FFFF  ""     
##  6 15669    6                Clifton Ro… ""            3 https:… 008080  ""     
##  7 15670    8                North Drui… ""            3 https:… 008000  ""     
##  8 15671    9                Boulevard … ""            3 https:… 00FFFF  ""     
##  9 15672    12               Howell Mil… ""            3 https:… FF00FF  ""     
## 10 15673    14               14th Stree… ""            3 https:… 00FF00  ""     
## # … with 108 more rows, and abbreviated variable names ¹​route_long_name,
## #   ²​route_desc, ³​route_type, ⁴​route_url, ⁵​route_color, ⁶​route_text_color
## # ℹ Use `print(n = ...)` to see more rows
## 
## $shapes
## Simple feature collection with 394 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -84.67085 ymin: 33.4323 xmax: -84.08274 ymax: 34.10663
## Geodetic CRS:  WGS 84
## First 10 features:
##    shape_id                       geometry
## 1    119402 LINESTRING (-84.55227 33.77...
## 2    119403 LINESTRING (-84.3694 33.821...
## 3    119404 LINESTRING (-84.39192 33.75...
## 4    119405 LINESTRING (-84.50486 33.73...
## 5    119406 LINESTRING (-84.21014 33.85...
## 6    119407 LINESTRING (-84.30635 33.88...
## 7    119408 LINESTRING (-84.33644 33.95...
## 8    119409 LINESTRING (-84.44857 33.65...
## 9    119410 LINESTRING (-84.35409 33.75...
## 10   119411 LINESTRING (-84.35409 33.75...
## 
## $stop_times
## # A tibble: 1,033,990 × 5
##    trip_id arrival_time departure_time stop_id stop_sequence
##    <chr>   <time>       <time>         <chr>           <int>
##  1 8044645 14:35:00     14:35:00       903320              1
##  2 8044645 14:36:02     14:36:02       903448              2
##  3 8044645 14:36:51     14:36:51       901144              3
##  4 8044645 14:37:43     14:37:43       904219              4
##  5 8044645 14:38:19     14:38:19       903850              5
##  6 8044645 14:38:39     14:38:39       903664              6
##  7 8044645 14:39:39     14:39:39       213268              7
##  8 8044645 14:40:12     14:40:12       213101              8
##  9 8044645 14:41:35     14:41:35       212926              9
## 10 8044645 14:41:53     14:41:53       211886             10
## # … with 1,033,980 more rows
## # ℹ Use `print(n = ...)` to see more rows

# Interactive mapping
tmap::tmap_mode('view')

## tmap mode set to interactive viewing

m1 <- tmap::tm_shape(atlsf$shapes) + tmap::tm_lines(alpha = 0.5)
m2 <- tmap::tm_shape(atlsf$stops) + tmap::tm_dots(id = 'stop_name', alpha = 0.5)

tmap::tmap_arrange(m1, m2, sync = T)

Calculate travel times from one station to all other stations

# Filtered stop times
stop_dist <- stop_group_distances(atlsf$stops, by='stop_name') %>% filter(dist_max > 200)

# Replace duplicate names of the stops by appending (1), (2), etc. so that 
# there will be no same stop_name
atl$stops <- atl$stops %>% 
  group_by(stop_name) %>% 
  mutate(stop_name = case_when(stop_name %in% stop_dist$stop_name ~ paste0(stop_name, " (", seq(1,n()), ")"),
                               TRUE ~ stop_name)) 

atl$transfers <- gtfs_transfer_table(atl, 
                                     d_limit = 200, 
                                     min_transfer_time = 120)

## ▶ Finding neighbouring services for each stop

## Loading required namespace: pbapply

## ✔ Found neighbouring services for each stop

## ▶ Expanding to include in-place transfers

## ✔ Expanded to include in-place transfers

am_stop_time <- filter_stop_times(gtfs_obj = atl, 
                                  extract_date = "2021-08-14",
                                  min_departure_time = 3600*7,
                                  max_arrival_time = 3600*10)
# travel_times
trvt <- travel_times(filtered_stop_times = am_stop_time,
                     stop_name = "MIDTOWN STATION",
                     return_coords = TRUE)

# Average time to get to other stations
print(str_c("On average, it takes ", 
            round(mean(trvt$travel_time)/60,1), 
            " minutes to travel from the Midtown Station to all other stations."))

## [1] "On average, it takes 50.1 minutes to travel from the Midtown Station to all other stations."

Visualize the routes and stops

GTFS from MARTA contains information about both rails and buses. The ‘routes’ table in atl has a column named ‘route_type,’ which contains integer values denoting types of transit service. To see what each integer means, you will need to look at the link to the Reference provided above (or see below for a screen capture from the Reference).

Reference from Google

‘routes’ table contains route_types, which is needed for seperating rail transit. ‘shapes’ table contains the shapefile needed for geooperations such as buffer, intersections, etc. Because ‘routes’ and ‘shapes’ table do not contain a common key, they need to be joined through ‘trips’ as intermediary table.

# Filter out rail transit
route_trip <- atl$trips %>% dplyr::left_join(atl$routes, by = "route_id")

unique_shape <- route_trip %>% 
  dplyr::group_by(route_id) %>% 
  dplyr::slice(1) %>% 
  dplyr::ungroup()

route.shape <- atlsf$shape %>% inner_join(unique_shape, by = "shape_id")

# Route type is not really intuitive - let's fix that
route.shape <- route.shape %>% 
  dplyr::mutate(route_type = dplyr::case_when(
    route_type == "0" ~ 'Tram, Streetcar',
    route_type == "1" ~ 'Subway, Metro',
    route_type == "2" ~ 'Rail',
    route_type == "3" ~ 'Bus'
  ))

pal <- leaflet::colorFactor(c("red", "orange", "pink"), domain = route.shape$route_type)

route.shape %>% 
  leaflet::leaflet(data = .) %>% 
    leaflet::addProviderTiles(providers$CartoDB.DarkMatter) %>% 
    leaflet::addPolylines(color = ~pal(route_type), 
                 weight = 3,
                 opacity = 0.9,
                 popup = paste0("Route type: ", route.shape$route_type))

Until Version 1.0, the sf package used ‘equirectangular projection’ for many operations that involve, for example, intersects, intersection, union, nearest_point, st_join, etc. It did so even when the given data is in geographical coordinates.

From Version 1.0, sf package uses S2 spherical geometry for spatial operations, based on S2 geometry library written by Google. For more information on what all of these means, you can read this post.

Calculate area of service by SES

GTFS from MARTA contains information about both rails and buses. The ‘routes’ table in atl has a column named ‘route_type,’ which contains integer values denoting types of transit service. To see what each integer means, you will need to look at the link to the Reference provided above (or see below for a screen capture from the Reference).

Reference from Google

‘routes’ table contains route_types, which is needed for seperating rail transit. ‘shapes’ table contains the shapefile needed for geooperations such as buffer, intersections, etc. Because ‘routes’ and ‘shapes’ table do not contain a common key, they need to be joined through ‘trips’ as intermediary table.

# Filter out rail transit
route_trip <- atl$trips %>% dplyr::left_join(atl$routes, by = "route_id")

unique_shape <- route_trip %>% 
  dplyr::group_by(route_id) %>% 
  dplyr::slice(1) %>% 
  dplyr::ungroup()

route.shape <- atlsf$shape %>% inner_join(unique_shape, by = "shape_id")

# Route type is not really intuitive - let's fix that
route.shape <- route.shape %>% 
  dplyr::mutate(route_type = dplyr::case_when(
    route_type == "0" ~ 'Tram, Streetcar',
    route_type == "1" ~ 'Subway, Metro',
    route_type == "2" ~ 'Rail',
    route_type == "3" ~ 'Bus'
  ))

pal <- leaflet::colorFactor(c("red", "orange", "pink"), domain = route.shape$route_type)

route.shape %>% 
  leaflet::leaflet(data = .) %>% 
    leaflet::addProviderTiles(providers$CartoDB.DarkMatter) %>% 
    leaflet::addPolylines(color = ~pal(route_type), 
                 weight = 3,
                 opacity = 0.9,
                 popup = paste0("Route type: ", route.shape$route_type))

Until Version 1.0, the sf package used ‘equirectangular projection’ for many operations that involve, for example, intersects, intersection, union, nearest_point, st_join, etc. It did so even when the given data is in geographical coordinates.

From Version 1.0, sf package uses S2 spherical geometry for spatial operations, based on S2 geometry library written by Google. For more information on what all of these means, you can read this post.

# If we do not convert CRS to a projected one, the buffer may generate a ragged boundaries. 
# https://r-spatial.github.io/sf/articles/sf7.html#buffers-1
# https://r-spatial.org/r/2020/06/17/s2.html#sf-10-goodbye-flat-earth-welcome-s2-spherical-geometry
sf::sf_use_s2(TRUE)
MARTA_buffer <- route.shape %>% 
  sf::st_transform(crs = 26967) %>% 
  sf::st_buffer(dist = units::set_units(400, "m")) 

MARTA_buffer_group <- MARTA_buffer %>% 
  group_by(route_type) %>% 
  summarise()

pal <- colorFactor(palette = c("red", "yellow", "blue"), domain = MARTA_buffer_group$route_type)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = MARTA_buffer_group %>% st_transform(crs = 4326), col = ~pal(route_type),
              popup = MARTA_buffer_group$route_type) 

Overlaying GTFS with Census

acs2020c <- acs2020 %>% 
  select(GEOID,
         hhinc = hhincE,
         r_tot = r_totE,
         r_wh = r_whE,
         r_bl = r_blE,
         tot_hh = tot_hhE,
         own_novhc = own_novhcE,
         rent_novhc = rent_novhcE) %>% 
  mutate(pct_wh = r_wh / r_tot,
         pct_bl = r_bl / r_tot,
         pct_novhc = (own_novhc + rent_novhc)/tot_hh) %>% 
  mutate(area1 = unclass(st_area(.))) %>% 
  st_transform(26967) %>% 
  mutate(area2 = unclass(st_area(.))) %>% 
  st_transform(crs = 4326) %>% 
  mutate(ln_pop_den = log((r_tot / (area1/1000^2)) + 1)) %>% 
  filter(!is.na(hhinc), !is.na(r_tot), !is.na(own_novhc))

pal1 <- colorNumeric(palette = "YlOrRd", domain = acs2020c$hhinc)
pal2 <- colorFactor("Spectral", domain = MARTA_buffer_group$route_type)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(group = "ACS",
              data = acs2020c, 
              color = "grey", 
              fillColor = ~pal1(hhinc), 
              fillOpacity = 0.5,
              weight = 1, 
              popup = leafpop::popupTable(round(st_drop_geometry(acs2020c[,c("hhinc", "pct_wh", "pct_bl")]),2))) %>% 
  addPolygons(group = "MARTA",
              data = MARTA_buffer_group %>% 
                st_transform(crs = st_crs(acs2020c)), 
              color = ~pal2(route_type),
              weight = 1,
              opacity = 0.9) %>% 
  addLayersControl(
    overlayGroups = c("ACS", "MARTA"),
    options = layersControlOptions(collapsed = FALSE)
  )

# Intersect buffer with tract 
buffer_per_tract1 <- acs2020c %>% 
  st_transform(crs = 26967) %>% 
  st_intersection(MARTA_buffer_group) %>% 
  filter(route_type == "Bus") %>% 
  # After Intersection
  mutate(subarea = unclass(st_area(.)),
         pct_served1 = subarea/area1,
         pct_served2 = subarea/area2) %>% 
  st_transform(crs = 4326)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

buffer_per_tract2 <- acs2020c %>% 
  left_join(buffer_per_tract1 %>% select(GEOID, pct_served1) %>% st_set_geometry(NULL),
            by = "GEOID") %>% 
  # There are many NAs in the pct_served1 column
  mutate(pct_served1 = case_when(is.na(pct_served1) ~ 0,
                                 TRUE ~ pct_served1))

pal <- colorNumeric(palette = "YlOrRd", domain = buffer_per_tract2$pct_served1)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = buffer_per_tract2, 
              fillColor = ~pal(pct_served1), fillOpacity = 0.9, 
              color = "white", opacity = 0.2, 
              weight = 1) %>% 
  addLegend("bottomright",
            pal = pal,
            values = buffer_per_tract2$pct_served1,
            title = "% Area within 400m from transit line") 

Be careful when calculating area. tm_shape( buffer_per_tract %>% mutate(diff = (area2 - area1)/area2)) + tm_polygons(col = 'diff', style = 'quantile')

# Is there any correlation between SES and service area?
var_name <- list(
  'hhinc'="Annual Household Income",
  'pct_bl'="% Black population",
  'pct_novhc'="% Household without vehicle",
  'pct_wh'="% White population",
  'r_tot'="Total population",
  'ln_pop_den' = "Population density"
)

var_labeller <- function(variable, value){
  return(var_name[value])
}

na_index <- buffer_per_tract2 %>% 
  select(hhinc, r_tot, ln_pop_den, starts_with('pct')) %>% 
  is.na(.) %>% 
  apply(., 1, function(x) sum(x)==0)

buffer_per_tract2 %>% 
  filter(na_index) %>% 
  pivot_longer(cols = c('hhinc', 'pct_wh', 'pct_bl', 'pct_novhc', "r_tot", "ln_pop_den"), names_to = "variable", values_to = "value") %>% 
  ggplot() +
  geom_point(aes(x = pct_served1, y = value), alpha = 0.4) +
  facet_wrap(~variable, scales = "free_y", labeller = var_labeller) +
  labs(x = "% Area with 400m beffer from transit line") +
  theme_bw()

## Warning: The labeller API has been updated. Labellers taking `variable` and
## `value` arguments are now deprecated. See labellers documentation.

test_var <- c("hhinc", "pct_bl", "pct_wh", "pct_novhc", "r_tot", "ln_pop_den")

map_chr(test_var, function(x){
  cor.object <- cor.test(buffer_per_tract2[['pct_served1']], 
                         buffer_per_tract2[[x]])
  r <- round(cor.object$estimate,3)
  t <- round(cor.object$statistic,3)
  p <- round(cor.object$p.value,3)
  return(paste0("for ",x ,", r=", r, " (t=", t, ", p=", p, ")"))
})

## [1] "for hhinc, r=-0.403 (t=-10.691, p=0)"   
## [2] "for pct_bl, r=0.089 (t=2.16, p=0.031)"  
## [3] "for pct_wh, r=-0.106 (t=-2.575, p=0.01)"
## [4] "for pct_novhc, r=0.52 (t=14.761, p=0)"  
## [5] "for r_tot, r=-0.307 (t=-7.819, p=0)"    
## [6] "for ln_pop_den, r=0.54 (t=15.553, p=0)"

Reflecting the frequency differences

Since some stops have more frequent services than others, it is useful to consider the service frequency. Here, we are doing to focus on the morning commute time (7AM to 10 AM) and count the number of departures at each stops.

service_ids <- atl$.$dates_services %>% pull(service_id)
stop_freq <- get_stop_frequency(atl, start_time = 3600*7, end_time = 10*3600, service_ids = service_ids, by_route = T) 
stop_freq_sf <- atlsf$stops %>% 
  left_join(stop_freq, by="stop_id") %>% 
  filter(!is.na(n_departures))

freq_pal <- colorNumeric("Reds", stop_freq_sf$n_departures)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addCircles(data = stop_freq_sf,
             fillColor = freq_pal(stop_freq_sf$n_departures),
             fillOpacity = 0.8,
             weight = 10,
             stroke = FALSE,
             popup = str_c(stop_freq_sf$stop_name, ": ", stop_freq_sf$n_departures)) %>% 
  addLegend(position = "bottomright",
            pal = freq_pal,
            values = stop_freq_sf$n_departures,
            title = "Count of Departures <br> between 7AM and 10AM")

stop_tract <- acs2020c %>% 
  st_transform(crs = 4326) %>% 
  st_join(stop_freq_sf) %>% 
  group_by(GEOID) %>% 
  summarise(n_departures = sum(n_departures),
            hhinc = mean(hhinc),
            r_tot = mean(r_tot),
            pct_wh = mean(pct_wh),
            pct_bl = mean(pct_bl),
            pct_novhc = mean(pct_novhc),
            r_tot = mean(r_tot),
            ln_pop_den = mean(ln_pop_den))

knitr::kable(stop_tract)

z <- function(vector){ round((vector - mean(vector))/sd(vector), 3) }

pass_per_stop_pal <- colorQuantile(palette = "Reds", domain = stop_tract$n_departures)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = stop_tract, 
              fillColor = pass_per_stop_pal(stop_tract$n_departures),
              fillOpacity = 0.8,
              color = 'grey',
              weight = 1, 
              popup = paste0("ID: ", stop_tract$GEOID, ", Value: ", stop_tract$n_departures)) %>% 
  addLegend(position = "bottomright",
            values = stop_tract$n_departures,
            pal = pass_per_stop_pal)

message(str_c("there are ", sum(is.na(stop_tract$n_departures)), " NAs in n_departures"))

## there are 84 NAs in n_departures

plotly::ggplotly(departure_count)